Two-State Migration of DNA in a Structured MicroChannel 
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DNA migration in topologically structured microchannels with periodic cavities is investigated 
experimentally and with Brownian dynamics simulations of a simple bead-spring model. The results 
are in very good agreement with one another. In particular, the experimentally observed migration 
order of A- and T2-DNA molecules is reproduced by the simulations. The simulation data indicate 
that the mobility may depend on the chain length in a nonmonotonic way at high electric fields. This 
is found to be the signature of a nonequilibrium phase transition between two different migration 
states, a slow one and a fast one, which can also be observed experimentally under appropriate 
conditions. 



I. INTRODUCTION 

DNA electrophoresis is one of the main techniques to 
separate DNA molecules by size Q . Since the mobility of 
DNA molecules does not depend on the chain length in 
free solution, electrophoresis is usually performed in gels. 
Given the progress of microtechnology, much effort has 
been spent on the integration of DNA electrophoresis into 
microfluidic devices (lab-on-a-chip devices) yl and sepa- 
ration strategies which rely on well-defined microscopic 
structures have been explored. Various such gel-free de- 
vices exploiting both DC [l|l|l|E[ll8| and AC |l[l3 
fields have been proposed. The mobility of the DNA in 
these microstructures is determined by a subtle interplay 
of the characteristic dimensions of the microstructure and 
the molecular size. Microstructured systems allow for 
fast separation of kbp size DNA fragments in time scales 
less than one minute [3, llfl Ull ■ 

In the present work, we focus on a geometry which has 
been investigated by Han and Craighead et al. 0, [l3i H^ 
and later by us (Duong et al. jg). The DNA is driven 
through a periodic sequence of cavities and constrictions 
(Fig. ^1 with an electrical DC field. The separation ex- 
periments of Han and Craighead motivated several com- 
puter simulation studies [ij, [l^ [T^ , which reproduced 
the results and lead to an improved understanding of the 
DNA separation mechanisms in these channels. Origi- 
nally [3, UM I the width of the constriction was chosen of 
the order of 0.1/im (corresponding to a few persistence 
lengths of DNA molecules), and entropic trapping at the 
entrances of the constrictions was proposed to be the 
main mechanism leading to size dependent migration of 
DNA in these channels. Recently, an additional effect 
that relies on diffusion has been revealed in such geome- 
tries |l6| . This second mechanism does not premise ex- 
tremely narrow constrictions. Therefore, one should also 
achieve DNA separation with wider channels. 

Based on that prediction, we have investigated experi- 
mentally the size-dependent electrophoresis in structured 
microchannels (Fig. ^| with constriction widths corre- 
sponding roughly to the radius of gyration of the tested 
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FIG. 1: (a) SEM image of an experimental device presented 
by Duong et al. ^, and (b) schematic drawing of the channels. 
As all size parameters are equal, only one parameter H is 
necessary to describe the structure. The coordinate system 
used in the simulations is given on the right. The dotted 
lines represent the electric field lines. Also shown in (b) is 
a configuration snapshot from a simulation of a chain with 
A^ = 320 monomers. 



DNA molecules |^ [Tl| . We found that separation by size 
is possible ill!|. Unexpectedly, however, the relation be- 
tween mobility and chain length turned out to be inverse 
to that observed by Han and Craighead: the short A- 
DNA (48 kbp) migrated faster than the longer T2-DNA 
(164 kbp)! 

In the present paper, we explore the reasons for this 
unforeseen behavior by Brownian dynamics simulations 
of a simple bead-spring model. Models of this kind have 
been used successfully in the past to study the migration 
of DNA in various geometries 11 IH El El llfl . Our 
simulations reveal a surprisingly complex phenomenol- 
ogy. If the model parameters are adjusted to those of 
the experiment, we reproduce the experimental findings. 



With other model parameters (smaller structures, lower 
electric field), we recover the originally expected increase 
of mobility with increasing chain length. As we show 
in this paper, the reversal of this trend turns out to be 
the signature of a nonequilibrium first order phase transi- 
tion, which emerges at even higher fields or in even larger 
structures. However unfavourable for separation in spec- 
ified geometries, this prediction implies the existence of 
two migration states, which could indeed be observed in 
migration experiments. 

The paper is organized as follows: the simulation 
model and the simulation technique are described in 
Sec.m and the experimental setup in Sec. lIIII In Sec. lIVI 
we briefly discuss the simulations corresponding to the 
earlier experimental study by Duong et al. 8]. The re- 
sults on the two coexisting migration states arc presented 
in Sec. and the possible transition mechanisms are 
discussed in Sec. IVII We summarize and conclude in 
Sec.lVnl 



II. THE SIMULATION MODEL 

A. Definition of the model 

We describe the DNA on a coarse-grained level by a 
bead-spring model. Several base pairs of DNA are rep- 
resented by one bead or "monomer" . Neighboring beads 
are connected by harmonic springs 



V.Ar) = 2^2, 



(1) 



where r is the distance between the two neighboring 
monomers. All beads interact with a purely repulsive 
Weeks- Chandler- Andersen (WCA) potential 



Kai^W/fesT: 



^Ufy'-ifr + \ ■■ (^)<2V6, 







otherwise, 



(2) 

where r is the distance between two monomers and a — 1 
sets the range of the interaction. 

The walls are assumed to be soft and purely repulsive. 
The potential describing the interaction with the walls is 
the same as that acting between the monomers (Eq. El 
with r as the the minimum distance between a bead and 
the wall). In the corners, the forces of the adjacent walls 
are summed up. To model the finite depth of the device, 
we also introduce walls in y-direction. 

Each bead carries a charge q and is subject to an elec- 
trical force /oi = qE — —qV^, but charges do not interact 
with one another in our model. The field is calculated 
by solving the Laplace equation (A<I> = 0) in the chan- 
nel. We use von Neumann boundary conditions at the 
walls {n ■ V(/> = 0, with ft the surface normal) and im- 
pose a constant potential difference between the inlets of 
the device. A numerical solution $(r) is obtained from 
the software program "Matlab" (Mathworks, US), using 
a finite element solver. The grid values of the derivatives 



are interpolated bilinearly. The field lines are shown in 
Fig. n Throughout the rest of this paper, the field de- 
noted by E refers to the total potential difference divided 
by the total channel length, and can thus be regarded as 
the characteristic strength of the electric field. 

The motion of the chain is described by Langevin dy- 
namics, i.e. the solvent surrounding the chain is replaced 
by a Brownian force ffi and an effective friction coefficient 
C = 1, which fulfill the conditions |2J 

<m> = (3) 

< V^At)VJAt') > = 2C keT 5,j5^0 S{t - t'), (4) 

with monomer indices i, j ^ 1 . . . N, cartesian directions 
a, /3 G {x, y, z} and t, t' two given times. The equations 
of motion are thus 



n = Vi 

mvi ^ fi - Qvi +rfi, 



(5) 
(6) 



where {fi} and {vi} represent the locations and the ve- 
locities of the monomers, respectively. The vectors {fi} 
sum up forces acting on a monomer (harmonic spring, 
repulsive monomer interaction, wall interaction, external 
electric field). Here, each bead carries a mass m, but 
simulations with vanishing mass were also carried out for 
comparison (see Sec. lVI|l . With m = 0, the equations of 
motion reduce to 



(n ^ fi + rfi 



(7) 



The natural units of our simulation are defined in terms 
of the bead size a, the friction coefficient (, the bead 
charge |g|, and the temperature T. Based on these quan- 
tities, the energy unit is k^T, the length unit is a, the 
time unit is to = (^a"^ /k-^T, and the electric field unit is 
Eo = k^T/a\q\. 

The dynamical equations (Eq. El and Eq. ^ are inte- 
grated with a Verlet algorithm using a timestep A, = 
lO^^to. In the case of the relaxational dynamics in 
Eq. 13 we use an Euler algorithm with the time step 
At = 10~^to- The stochastic noise f^i was implemented 
by pic king random numbers at every time step. As shown 
in Ref. |23, the random numbers do not have to be Gaus- 
sian distributed, as long as they fulfill Eq. (31 and Eq. 01 
Those used in our simulation were evenly distributed in- 
side the unit sphere. Unless stated otherwise, the run 
lengths were 4 — 20 • 10* At with equilibration times of 
2 • 10^ At. 

The spring constant in Eq. ^ was chosen very high, 
k = lOQk^T/a^ . This ensures that no chain crossing oc- 
curs (23. The equilibrium length of each bond is 0.847(7. 

The static properties of free chains of length N corre- 
spond to those of self-avoiding random walks (23i] with 
the persistence length l^ ~ 1.6cr and the radius of gyra- 
tion i?g « 0.5cr • iV. Here v = 0.588 is the Flory expo- 
nent. The dynamical properties are characterized by the 
diffusion constant D = k^T/NQ = N~^a^/tQ, the mo- 
bility /io — \q\/C — o'/Eoto, which is independent of the 



chain length A^ , and the decay time of the drift velocity 
tq = m/C, (Eq. EJ . The quantity Tq sets the time scale 
on which inertia effects are significant. In most simula- 
tions presented here, it was set to tq — to (i.e. m — C^o)j 
based on the assumption that the relevant time scales for 
the dynamical processes of interest here are much longer. 
In some cases, however, simulations with vanishing tq 
(vanishing mass) were also necessary (see Sec. IVI|I . On 
time scales larger than tq, the chains behave like stan- 
dard (self-avoiding) Rouse chains |22 with the Rouse re- 
laxation time Tr « 0.045 to ■ N^+'^''. 



B. Correspondence bet^veen the model and the 
experiment 

In the experiments, the mobility of A-DNA (48 kbp) 
and T2-DNA (164 kbp) in microchannels with geometries 
such as sketched in Fig. ^ was studied. Details are given 
in Ref. la and Sec. IIIII We will now discuss how the 
parameters and units of our model can be related to those 
of the experiment. 

The energy unit is simply given by the temperature, 
keT = SOOfce K « 0.026eV, at which the experiments 
were carried out. Next we adjust the value of the length 
unit. In our earlier work on entropic traps |la |. we had 
matched the persistence lengths of the model chains with 
that of DNA. In the present case, however, the per- 
sistence length is not an experimentally relevant length 
scale, because all dimensions of the microchannels are on 
the order of magnitude of the radius of gyration. Thus we 
can determine the length unit by comparing the charac- 
teristic length H in our model channel (see Fig.^ to the 
corresponding experimental value. Note that this implies 
that we can relate simulations with one channel geometry 
to all three experimental channels by scaling our results 
accordingly. 

More specifically, the simulations were carried out us- 
ing channels with H = 60cr. We illustrate our adaptation 
with the example of the 5/xm microchannels. It is carried 
out in four steps. First, we compare the channel height 
H = bfim = 60a. This gives the length unit la = 83nm. 
Second, comparing the radii of gyration of our model 
chains with that of A-DNA, we find that A-DNA is rep- 
resented by chains of A^ = 140 beads, or 1 bead = 340 
bp. Third, the time unit is adjusted by matching the 
diffusion constant D for A-DNA. Experimentally, Smith 
et. al. have reported D = 0.47 ±0.03/im^ /sec 24]. Com- 
paring this with the theoretical value, D = a'^ /Nto, we 
obtain to = 1.1 • 10~'*s. Hence, one second corresponds 
to 9.5 • lO'^to (or ^ 10^ time steps). Finally, the elec- 
tric field unit Eq is adjusted by adjusting the mobility in 
large unstructured microchannels, which is chain length 
independent and given by ^o = 1.84 • lO^^cm^/Vs [^. 
It is the sum of the free flow mobility /io and the elec- 
troosmotic mobility /i^^f — (2.9 ±0.6) • 10~''cm^/Vs in our 
channels (see below). This leads to IV/cm = 2.3T0~^£'o. 
These values as well as those obtained for channels with 



H — 3fj,m and H = 1.5/^m are summarized in table Q] 
Unless stated otherwise, all adaptations given in this pa- 
per refer to the 5//m channels, except in Sec. IIVI 

We set the depth of the device to 60cr. For compar- 
ison, we also performed simulations with infinite depth 
of the device and found no qualitative difference (data 
not shown). In our case, the depth has no significant 
influence on the dynamics. 

The model disregards a number of important physical 
effects. 

First, electrostatic interactions are neglected. The De- 
bye screening length of DNA in typical buffer solutions is 
about 2nm, which is comparable to the diameter of the 
polymer and much smaller than the persistence length of 
DNA (roughly 50nm). 

Second, PDMS exhibits silanol groups on its sur- 
face which dissociate under the experimental conditions. 
Thus the experiments were carried out in microchannels 
with negatively charged surfaces. This implies the gen- 
eration of cathodic electroosmotic flow, and the resulting 
mobility of DNA molecules is a sum of the electroosmotic 
and electrophoretic mobilities. The DNA molecules mi- 
grate to the anode, indicating that electrophoresis over- 
comes electroosmosis. In the simulation, we do not incor- 
porate any electroosmotic flow, or flow in general. How- 
ever, recent reports 25, 26] show that for sufficiently low 
Reynolds number, the flow outside the Debye layer at 
the walls is rotation free and proportional to the electric 
field with a fixed proportionality constant a. Let vo{r) 
denote the flow field at a given position r and /o,i the 
forces acting on monomer i, excluding the electric force 
fai,i = — 9- V$(ri). Using vo{r) = a- V<I'(r), the equation 
of motion (Eq. (BJ now becomes 



fa,l + fal.t ~ C{Vz + Vo{r)) + T}i 



(8) 



= /o,» - (g + Ca) • V$(r) + r/„ 

and for the overdamped dynamics (Eq. 01 one finds 

C'^»=£..-('? + Ca)-V$(rO+77,. (9) 

These equations are formally equivalent to Eq. El and 
Eq. [7| with a rescaled charge q. The Reynolds number 
for the electroosmotic fiow is easily found for our device. 





H = l.Spm 


H = ?,^m 


H — 5/irn 


1 bead 


48 bp 


150 bp 


340 bp 


A-DNA 


1000 beads 


330 beads 


140 beads 


T2-DNA 


3500 beads 


1200 beads 


490 beads 


1/im 


40cr 


20cr 


12a 


Is 


7.5 ■ lO'-'io 


6.2 • 10*to 


9.5 ■ lO^fo 


IV/cm 


9.8- IQ-'^Eo 


5.9 • lQ-*Eo 


2.3 • 10-^So 


£^0 


lOkV/cm 


1.7kV/cm 


430V/cm 



TABLE I; Adaptation of the simulation units to various chan- 
nel sizes. All values have been rounded to two digits. 



Using the viscosity of water, the electroosmotic mobil- 
ity given above and a characteristic length of 5/ini, the 
Reynolds number is roughly 10^'^ for electric fields of 
100 V/cm in our channels. This justifies our adaptation 
to the overall mobility given above. 

Third, hydrodynamic effects are not taken into ac- 
count. This approximation must be questioned. On 
the one hand, DNA is always surrounded by counter 
ions, which are dragged into the opposite direction of the 
DNA. Thus the DNA molecule experiences not only hy- 
drodynamic drag, but also an extra friction from the sol- 
vent molecules. In free solution, these two effects cancel 
each other [l| . This "hydrodynamic screening" accounts 
for the free-draining property of DNA. However, the to- 
tal cancellation fails if the DNA molecule is blocked by 
a geometric barrier [23|. In that case, the counter ions 
will not be immobilized, since counterions are still free 
to move. Furthermore, hydrodynamic interactions affect 
the diffusion constant D. In our model, it scales the 
chain length N like a Rouse chain {D oc 1/N). Including 
hydrodynamic interactions, one would expect Zimm scal- 
ing {D oc 1/Rg (X 1/N'^). In experiments, the diffusion 
constant of DNA is found to scale as D ex 1/N°-^'^'^ JH]. 

Unfortunately, a full simulation which treats hydrody- 
namic interactions correctly, takes into account the dy- 
namical counter ion distribution, and includes electro- 
static interactions, is a formidable task and impossible 
with standard computer resources. The simplifications 
described above will affect the quantitative results, but 
presumably neither by orders of magnitude nor qualita- 
tively. 



taining 50 mM NaCl, 1 mM ethylenediaminetetraacetic 
acid (EDTA) and 2% (v/v) /3-Mercaptoethanol) were 
adjusted to a YOYO-1 (Molecular Probes, USA) in a 
basepair ratio of 1 : 7.5. A sensitive fluorescence video- 
microscopy setup is used to record DNA migration, as 
described in Ref. [3. Briefly, it consists of an inverted 
microscope (Axivert 100, Zeiss, GER) with a filterset 
for fluorescence observation (XF 100-3, Omega, USA), 
a lOOx oil immersion objective (Plan Neofluar NA 1.3, 
Zeiss, GER) and a sensitive CCD-camera (Imager 3LS, 
LaVision, GER). Data acquisition and automated data 
analysis is performed in DaVis 5.4.4 (LaVision) using 
cross-correlation analysis. Electric fields are applied with 
power supplies from FUG (MCN 14-2000, GER). 

The electroosmotic mobility was determined according 
to the current monitoring method [291] . Linear channels 
of 20 /im width and depth were filled with water immedi- 
ately after assembly, which was exchanged by the 10 mM 
sample buffer without DNA before the experiment. The 
buffer in the anode reservoir was then exchanged by 8 mM 
Tris buffer with otherwise identical composition as the 
lOniM buffer. A constant electric field of 400 V/cm was 
applied and the current was recorded until a stable lower 
current value was reached. The recorded current decay 
times T^„f were used to calculate the electroosmotic mo- 



bility according to: /icof — 



ET„af 



-- ' JJJ /(jot 

channel, and E the applied electric field 



with / the length of the 



IV. LOW ELECTRIC FIELD: STEADY-STATE 
MIGRATION 



III. EXPERIMENTAL SETUP 

The fabrication method used for the structured 
microchannels is based on soft lithography of 
poly(dimethylsiloxanc)(PDMS) and is described in 
detail in Rcf. 8. Briefly, a photoresist (SU-8 from 
Microresist, GER) coated Si- wafer (CrysTec, GER) is 
exposed through a chromium mask (DeltaMask, NL). 
After developing and hard-baking liquid PDMS (Sylgard 
184, Dow Corning, USA) is poured onto the master 
containing the inverted microstructure and baked for 
3 h at 75°C. Peeling off the PDMS slab, punching 
reservoir holes through the PDMS and covering it with 
a clean microscope glass slide results in the microchip. 
Fig. 2b shows a scanning electron micrograph (SEM) 
image of the topview of a microchannel with H=5/im. 
In comparison to the simulated geometry, the corners 
of the cavities are rounded and the lengths a, b, c, and 
d in Fig. ^ are not exactly equal. The exact values 
were a=3.6/j,m, b= 6.0/im, c=3.7/im and d=4.7/im. The 
depth of the microchannel was 2.8/im. Other channel 
geometries mentioned in this work are described in detail 
in Ref. S 

For fluorescence imaging, 6pMA- or T2-DNA (Fluka, 
Germany) in the same buffer (10 mM Tris at pII8.3 con- 



We begin by presenting simulations for relatively low 
electric fields E, which can be compared to the previously 
published experimental results of Duong et al. [g. In this 
section, we will adjust our model units to the channels 
with characteristic size H — 3/xm (cf. Sec. IIIB|l . where 
Duong et al. had found particular clear evidence of a 
size-dependent mobility. 

The electric field was varied from E = 0.0025i?o up to 
E — 0.04_Eo, which corresponds to experimental fields in 
the range of 4.2V/cm up to 67V/cm. The chain lengths 
range from A^ = 10 up to A^ = 2000 monomers, modeling 
DNA strands from 1.5 kbp up to 290 kbp. Fig. |21 shows 
typical trajectories obtained at i? = 0.005i?o. The trajec- 
tories are similar to those obtained for the entropic trap 
geometry [l^. For N — 10, the trajectory is completely 
dominated by diffusion. Long chains with A^ ~ 1000 
monomers are hardly affected by the constrictions - much 
less, in fact, than in the case of the entropic traps |16i] . 
This is because the width of the narrow regions is of the 
order of the radius of gyration; the chain is thus able to 
cross the constriction without uncoiling. The trapping 
mechanism suggested by Han et al. p should not apply 
here. The fact that some trapping nevertheless occurs for 
short chains underlines the importance of the effect re- 
ported by Streek et al.jla|. Indeed, Fig.Oshows that the 
mobility increases as a function of the chain length for 
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FIG. 2: Trajectories for chains of length (a) N = 10, and (b) 
N = 1000, at _E = 0.005-Eo. The horizontal dashed hnes mark 
the beginning of the narrow regions. The solid line in the 
middle shows the position of the center of mass, the dashed 
lines indicate the leading and trailing monomers. For N = \Q 
(a), these lines cannot be distinguished from each other. 
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FIG. 3: Mobility as a function of chain length for E = 
0.005 So (8.4 V/cm) and E = 0.04 So (67 V/cm). The ver- 
tical dashed lines indicate the length of A— and T2-DNA, 
with adaptation for the 3/im constrictions. Experimentally, 
Duong et. al. have reported fi/fio ~ 0.99 ±0.1 for A-DNA and 
m/Mo = 0.63 ± 0.03 for T2-DNA |||. 



with larger structures. Therefore, we have fabricated and 
studied microstructures with H = 5/xni, and the results 
will be presented in the next section. 



E = 0.005i?o. In the experimental channel, E = 0.005i?o 
corresponds to a real field of _E ^^ 8.4V/cm. In such 
low fields, Duong et al. 8] were unable to distinguish 
between different mobilities for different DNA sizes; the 
statistical error is large because diffusive motion becomes 
significant. We note that in simulations, we can track the 
DNA migration over many more cavities than in the ex- 
periments, therefore the relative error is smaller. 

At higher fields, Duong et al. Q reported the inverse 
effect - the mobility of T2-DNA is smaller than that of A- 
DNA. If we increase the electric field in the simulations 
accordingly, we find indeed that, at E ^ O.OAEq, the 
mobility reaches a maximum at A^ = 200 beads (29 kbp) 
and decreases for longer chains (Fig.|2J). In the parameter 
region corresponding to A-DNA and T2-DNA, smaller 
chains migrate faster than longer ones, but the difference 
in the mobility is not as pronounced as in experiment. 
Thus our simulations reproduce the behavior observed 
by Duong et al. 8] qualitatively. 

This result is gratifying, but it does not yet explain 
why the chain length dependence of the mobility is sud- 
denly reversed. In order to explore this question, we 
proceeded to investigate the migration in our structures 
under more extreme conditions. In the simulations, we 
can increase the field strength E. In the experiments, 
this can only be realized up to a field of approximately 
100 V/cm, where the cross correlation analysis [8| fails. 
However, table ^ shows that we can alternatively work 



V. HIGH ELECTRIC FIELD: TWO MIGRATION 
STATES 

We begin by discussing the simulation results. We 
started with studying the extreme case E = Eq. This 
corresponds to 430 V/cm in a S/im structure and is thus 
not accessible in experiments. Two sample trajectories 
for a chain of length A^ = 10 (3.4 kbp) and N = 400 (140 
kbp) monomers are given in Fig. ^ As can be seen from 
the insets in the A^ = 10 case, the chains retain some 
memory of their state in the previous cavity. If a chain 
passes a constriction without being trapped, it will avoid 
trapping at the next barrier as well. A similar effect had 
already been observed in the entropic traps [Ig. As in 
the simulations discussed in the previous section, chains 
may still get trapped in the corners of cavities, but at a 
lower rate due to the reduced diffusion time inside a sin- 
gle cavity. For N = 400, the migration is periodically fast 
and slow, indicating that the polymer penetrates the low 
field regions in the cavities. Over long simulation times, 
both the trajectories for A^ = 10 and A^ — 400 arc quite 
regular. 

The situation is qualitatively different for A^ — 200 
(70 kbp). The trajectory is smooth on short time scales, 
but occasionally, the speed of migration changes abruptly 
(Fig. 0). One observes two states of migration, one of 
which is fast and the other is slow. The lower panel 
of Fig. O shows that the two states are associated with 
two distinctly different z-positions of the center of mass: 
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FIG. 4: Trajectories of a chain with a) A'^ = 10 (3.4 kbp) and 
b) iV = 400 (140 kbp) monomers at £ = Bo (430 V/cm). The 
dashed lines inside the inset in the N = 400 case represent 
the onset of the narrow regions. 



the polymer migrates faster when it stays in the homoge- 
neous field in the upper part of the channel. From Fig.|Sl 
both states are long-lived during the simulation. Their 
life time can be made arbitrary long by using large chan- 
nels (cf. Table P) . We conclude that our microchannel 
system exhibits a nonequilibrium first order phase tran- 
sition between two different migration states. We find 
coexistence regions of the two states (Fig.|Sland|Sl). The 
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FIG. 5: Trajectory of a chain with iV = 200 (70 kbp) 
monomers at _E = _Eo, along the x direction (top) and the z 
direction (bottom). Note the two different migration speeds 
which can be related to the penetration depth into the wide 
cavities (the dashed vertical lines are just guides to the eye) . 
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FIG. 6: Snapshots obtained from (a) an experiment with T2- 
DNA and (b) simulation with 200 beads a.t E = Eo. The 
experimental snapshots show a series taken with 40 ms time 
step (fast state shown by arrows pointing down, slow state 
indicated by arrows pointing up) . Both snapshots show chains 
in the slow state, which penetrate the wide region, form a 
coil inside and stretch when passing into the narrow region. 
Chains in the fast state are also shown. Here the polymer 
remains permanently coiled in the homogeneous part of the 
field. 



nonmonotonic chain length dependence of the mobility 
aX E — 0.04£'o is a consequence of that transition. 

The experimental studies of A-DNA and T2-DNA mi- 
grating in large structures support our simulation results. 
In structures characterized by iJ w 5/im, at electric fields 
ranging up to 100 V/cm, both A- and T2-DNA exhibit 
two distinct states of migration, which are illustrated in 
Fig. EJi. In one state, the chain remains coiled all the 
time in the upper (homogeneous) part of the field, and 
travels at high speed. In the other, it penetrates deep 
into each cavity, forms a coil inside, and gets stretched 
again when moving into the narrow region. 

Fig. 13 shows monomer density histograms obtained 
from experiment and simulation. In both the exper- 
iments and in the simulations with chains of length 
N = 120 — 250 beads, we observe a low-density depletion 
zone between the upper and the lower part of the channel. 
Furthermore, Fig. [7|i, shows that short chains are more 
likely to migrate in the fast state and long chains are more 
likely to migrate in the slow state. For N = 120 — 250, 
the polymer migrates alternatingly in the slow and fast 
state, as has already been seen from the trajectory of 




a) 



Length of DNA [kbp] (5 |im adaptation) 
10 100 



FIG. 7: (a)Monomer density histograms obtained from sim- 
ulation for N = 50,100,120,140 (top, left to right) and 
N = 170, 200, 250 and 320 (bottom) at E ^ Eq. The density 
distribution p was cut off at p = 1/2 to dampen the peaks 
at the corners of the constrictions and to emphasize the de- 
pletion zone, (b) Monomer density histograms obtained from 
experiment at _E = 86V/cm (0.20-Eo) for A-DNA in 5/im con- 
strictions. Note the depletion zone between both states. 



iV = 200 in Fig. m 

Density histograms such as those shown in Fig. [7^ al- 
low to determine the population density of the two states. 
We choose the z-position of the center of mass as an in- 
dicator of a polymer's state (cf. also Fig. |SJ|. Since the 
transition does not take place instantaneously (Fig. O, 
we define two thresholds separated by a gap: a chain is 
taken to be in the slow state or in the fast state, if z < Sa 
or z > 20ct over at least 5 • 10^ time steps. After having 
assigned a state to large portions of the trajectory, we 
can determine the mobility in the slow and the fast state 
separately. 

Fig. |SJi shows the results for E — Eq and compares 
it with the overall mobility. The assignment of states 
works well for intermediate and long chains. For N — 100 
(34 kbp) , the polymer migrates almost completely in the 
fast state and for polymers with N > 320 (110 kbp), the 
polymer occupies the slow state only. Short chains diffuse 
very strongly in the z direction, no depletion zone occurs, 
and they often cross from the lower to the upper part of 
the channel. Therefore, the assignment criterion often 
fails, and only a few short fragments of the trajectory 
contribute to the calculation of the mobility in the fast 
state. This explains why the result does not coincide 
with the total mobility at low A^. 

The transition was not only observed e.t E — Eq, but 
also observed at smaller fields E. Fig. ^ and Fig. ^ 
gives the mobilities and the population densities for 
E = 0.25 £:o (llOV/cm) and E = OM Eq (17V/cm), 
as determined by the above criterium. At E = 0.25 Eq, 
the transition is not as pronounced as for E = Eq, but 
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FIG. 8: Mobility as a function of the chains length at E = 
Eo (a), E = 0.25 Eq (b) and E = 0.04 Eo (c). Note the 
decrease in mobility at A = 100 - 250 (34 kbp - 85 kbp) at 
E = Eo, which is also reflected by the population densities 
given in the lower panel. For short chains, diffusion dominates 
the migration, and most of the time, it is not possible to assign 
unambiguously a single migration state. This is reflected by 
the amount of "undefined" population. The error bars of the 
slow and fast mobilities show the statistical error, based on 
the spread of the values. Note that the transition fades away 
when the electric field is decreased. 



still present. At E = 0.04 _Eo, the transition has almost 
disappeared, and only a slight decrease in the overall mo- 
bility remains. This confirms our earlier assertion that 
the nonmonotonic chain length dependence of the mobil- 
ity a,t E = 0.04 £^0 is a signature of the phase transition. 
It also clarifies why long chains (T2-DNA) migrate slower 
than shorter chains (A-DNA) in our channels. They tend 
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to penetrate deeper into the low field regions, which slows 
them down significantly. 

The quantitative comparison between simulations and 
experiments raises a question. According to the simu- 
lations, coexistence of two states should be observed at 
chain lengths below N — 250 beads, which corresponds 
to 64 kbp in 5/im microchannels. Thus we would expect 
two mobilities for A-DNA, but only one for T2-DNA. Ex- 
perimentally, we find that T2-DNA exhibits two states of 
migration as well. This is in contrast to the prediction 
of our simulation. However, a closer inspection reveals 
that even chains of A'^ = 500 beads (170 kbp, roughly 
T2-DNA) remain in a fast state for a long time (up to 
8 • 10* time steps, data not shown), if they are prepared 
accordingly. These initial fast states turned out to be so 
stable, that we were unable to analyze their life time in 
detail. The experimental situation hardly corresponds to 
an "equilibrated" late-time limit. It is conceivable that 
after being introduced into the microfluidic channel, a 
sizable fraction of molecules happens to be prepared in 
a state that is very weakly populated in the long-time 
average. 



VI. MIGRATION STATE TRANSITION 
MECHANISMS 

To understand the transition and identify the factors 
which stabilize the two migration states, we have investi- 
gated the transition processes between the two states in 
more detail. 



A. Transition from the fast to the slow state 

We begin with discussing the transition from the fast 
to the slow state. A series of snapshots of a chain with 
N = 500 beads in the process of crossing from the fast 
to the slow state is shown in Fig. IHl The chain travels 
through roughly 35 cavities during the transition process, 
which we can divide into five stages: 

1. At the beginning, the polymer is in a coiled con- 
formation in the homogeneous field region. The 
radius of gyration of a chain of A^ = 500 beads in 
free solution is R^ « 18a. The size of the homoge- 
neous field region is limited, thus larger polymers 
are deformed (Fig. |5Jl) . 

2. At some stage, one polymer loop or end happens to 
penetrate into a region with a lower field. As that 
polymer part experiences a lower pulling force, it 
slows down, while the force acting on the rest of 
the polymer is unchanged. This leads to uncoiling 
and stretching of the polymer. 

3. More and more monomers are slowed down. Since 
the field in the homogeneous part pulls horizontally 
in the x-direction, and the delayed monomers are 



located in the cavities at low z values, the total 
force acting on the chain has a net downward com- 
ponent (Fig. El b). After a while, the polymer is 
completely stretched (Fig. I^c). 

4. For the highly stretched chain, (Fig.^c), the prob- 
ability that the foremost monomer gets trapped in 
a cavity is very high. Immediately afterwards, the 
polymer collapses in that cavity. (Fig.^d). 

5. At the end, the polymer migrates in the slow state. 
For long polymers, the probability to leave the slow 
state is very low (Fig. Eld). 



B. Transition from the slow to the fast state 

Unfortunately, the configuration snapshots of the tran- 
sition from the slow to the fast state cannot be inter- 
preted as easily as those in Fig. O and we have not yet 
been able to distil the mechanisms which stabilize the 
fast state in a satisfactory way. 

One effect which presumably favors the transtion is 
the inertia of the molecule at the corners of the con- 
striction. Monomers migrating up at high speed against 
the wall continue to drift upwards slightly as they en- 
ter the constriction (data not shown). This effect in- 
volves dynamical processes on time scales smaller than 
To, the characteristic decay time for the drift velocity of 
the polymer. In our model, tq is given by tq — m/(^, 
which corresponds to approximately 10~^s for beads of 
mass m — C,tQ. For DNA in free solution, the drift veloc- 
ity decay time is [S^ tq « 10~^ . . . 10~^^s, which is much 
smaller. Thus this particular drift effect is an artifact of 
the simulation model. To verify that this mechanism is 
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FIG. 9: Transition of chains with N = 500 beads (170 kbp) 
from the fast to the slow state. The snapshots are shifted 
by the number of cavities given on the right. The transition 
itself takes about 35 cavities. In b) a force parallelogram is 
sketched schematically. The electric force F^ and the entropic 
spring force F^ on the leading monomers add up to a total 
force which has a net downward component. 





.l,f~«->» 






. ^ 




e 


c 


' 


^^__^ 






^ 



not responsible for the transition observed in the simu- 
lations, we modified the model by changing its dynamics 
to that of an overdamped system [211 (Eq. CJ . We could 
only simulate short chains, and were unable to perform a 
thorough analysis due to the reduced time step (Eq. [TJ. 
The transition still occurs in both directions (data not 
shown). This proves that the fast state is stabilized by 
another mechanism. In fact, the only noticeable effect of 
inertia was to suppress the transition for short chains. 

Another conceivable mechanism which might favor and 
stabilize the fast phase is based on diffusion. While the 
polymer migrates along the wall in the narrow region and 
around the corners, it can only diffuse upwards. This ef- 
fect is enhanced by the tendency of the polymer, which 
has been forced into a stretched conformation in the in- 
homogeneous parts of the field, to contract back to a coil. 
The characteristic time scale for contraction is presum- 
ably related to the Rouse time and scales like A^i+^'' |2ll| . 
which would explain why the slow state dominates at 
large chain length. We have not yet been able to support 
this explanation quantitatively. 

The understanding of the fast phase remains an open 
problem. Systematic studies are difhcult, because the 
opposite transition is triggered by diffusion as well. 



C. Critical chain length 

Fig. |51 suggests that the transition can only take place 
if the electric field exceeds a critical value; the inhomoge- 
neous field has to uncoil the chain and thus work against 
the entropic spring force of the polymers. The mean elon- 
gation (r) of stretched self-avoiding polymers subject to 
a force /, which pulls on the end monomers, is 



(r) oc Na {fa/k^Tf'^ 
for strongly stretched chains, and 

(r) ex a'N^-'/{k^T)f 



(10) 



(11) 



for weakly stretched chains, where a ~ 0.5a is the propor- 
tionality constant in i?^ = aTV. In our case, the force of 
course acts on all monomers. Nevertheless, we shall use 
Eg.llOland Eq.^Jto roughly estimate the maximal effect 
of the inhomogeneous field on a chain, with f (x E. 

The question is: how much does the field have to 
stretch a chain in order to stabilize a distinct slow state? 
Let us consider the critical chain length A^„ at which both 
states are equally populated. Assuming that it is suffi- 
cient to stretch the chains by an amount proportional to 
their own size, one would expect NE^^^ to be a constant 
at iV = iV, (according to both Eq. [TUland Eq. [TTl)- Al- 
ternatively, it might be necessary to stretch the chain to 
a fixed length, which is determined by the size of the de- 
vice, in particular the constriction length H (see Fig.^). 
In that case, according to Eq. ^| NE'^^^ should be con- 
stant. This is tested in table HU The quantity N, ■ E^^^ 
is obviously not constant for different electric fields E. 



E/Eo 


A^c 


iVe ■ {E/Eof^'' 


TVo ■ {E/Eof/'' 


0.08 


290 ± 25 


4.3 ±0.4 


55 ±5 


0.15 


174 ± 20 


7.4 ±0.9 


50 ±6 


0.25 


124 ± 15 


12.3 ±1.5 


50 ±6 


0.35 


120 ± 15 


20.9 ± 2.6 


60 ±8 


0.50 


106 ± 10 


33.4 ±3.2 


67 ±6 


0.70 


125 ± 10 


71.7 ±5.7 


103 ±8 


1.00 


200 ±5 


200.0 ±5.0 


200 ±5 



TABLE II: Crossover lengths A'^ for various electric fields. 
Errors of the crossover lengths are estimates of the population 
density histograms. 



However, N^ ■ E'^/^ is almost constant for E < 0.50i?o- 
For larger values of E this rule breaks down. This might 
be related to the inertia effects discussed in Sec. IVIBI 
Simulations with overdamped dynamics at i? = -Eq in- 
deed suggest that, as m — > 0, chains as small as iV = 50 
already exhibit two distinct migration speeds (data not 
shown). This would yield N^E'^/'^ = SOSq/^ which is 
consistent with results at lower fields. Unfortunately, the 
statistical quality of our data is not sufficient to deter- 
mine N^ rigorously in the overdamped case. 

Nevertheless, the sum of our findings suggests that the 
critical length is determined by the geometry of the de- 
vice rather than by the size of the chain. This could 
explain why the transition disappears at small fields: for 
E = 0.04i?o, the predicted critical length is N„ « 500, 
but free chains of that length have an end-to-end dis- 
tance of 46cr, which is comparable to the channel height 
H . Hence a fast state, in which the chain is essentially 
unperturbed, cannot exist. 

This result suggests strategies for the design of mi- 
crostructures which do or do not exhibit the two-state 
behavior; we expect it can be suppressed by choosing 
the constriction length d in Fig. ^p small, and it can be 
promoted by making d larger. 



VII. CONCLUSION 

We have presented Brownian dynamics simulations 
and experimental studies of DNA migration in struc- 
tured microchannels. Our simulations reproduce earlier 
experimental results presented by Duong et al. ^] and, 
in particular, explain the experimentally observed migra- 
tion order of A- and T2-DNA in migration studies and 
separation experiments S]: in channels with geometries 
as shown in Fig. ^ and for moderate field values, the 
shorter A-DNA molecules migrate faster than longer T2- 
DNA molecules. 

This behavior, which is opposite to that expected 
at weaker fields, is a signature of a high-field non- 
equilibrium phase transition. At very high fields (or in 
much larger structures), we found that chains can mi- 
grate with two distinctly different speeds, assuming dif- 
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ferent states of migration. This was observed both ex- 
perimentally and in simulations. We have discussed the 
factors which stabilize the two states and can thus pro- 
pose strategies for the design of channels which do or do 
not exhibit the transition, and which could be suitable 
for separation experiments. 

The comparison of simulations and experiments 
demonstrates that our simple DNA model, which dis- 
regards electrostatic interactions, hydrodynamic interac- 
tions, and irrotational fluid flow, nevertheless reproduce 
the experimental results of DNA migration in our mi- 
crochanncls. However, the interplay of DNA motion and 
buffer flow is non-trivial and may be neglected for irro- 
tational EOF outside Debye layers only. As mentioned 
in Sec. IIIBI the experimental electroosmotic mobility in 
straight channels is of the same order of magnitude as 
the DNA mobility. In structured microchannels, where 
DNA molecules migrate at either high fields or along the 
Debye layer, electroosmosis should give rise to interest- 
ing nontrivial flow patterns which certainly influence the 
DNA migration. Moreover, the electrostatic and the hy- 
drodynamic interactions are not entirely screened in mi- 
crochannels, which leads to additional effects [23|. In 
order to investigate such phenomena, systematic experi- 
mental studies are necessary, and efficient new simulation 
methods need to be developed. These shall be explored 
in the future. 

From a practical point of view, our results have two 
important implications: 

On the one hand, they demonstrate that the physics 



of DNA electrophoresis in microchannels is surprisingly 
complex. When designing new microfluidic devices, one 
must be aware of the possibility of nonmonotonic or even 
bistable behavior, and a thorough theoretical analysis is 
advisable. 

On the other hand, the observed two-state behavior 
can presumably be exploited. Our findings indicate that 
the migration velocities of the DNA fragments differ by a 
large factor (almost 2 in Fig.lHJ for N > N, and N < N^. 
Since the chain length corresponding to N^ depends on 
the channel geometry (cf. Table P), a device with a gra- 
dient in the constriction length could in principle sequen- 
tially separate molecules very efficiently. However, as the 
switching time between the two states is very long, meth- 
ods need to be developed to accelerate the transition. We 
are currently exploring corresponding strategies experi- 
mentally as well as theoretically. 
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